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ABSTRACT 

Various radio observations have shown that the hot atmospheres of galaxy clusters are mag- 
netized. However, our understanding of the origin of these magnetic fields, their implications 
on structure formation and their interplay with the dynamics of the cluster atmosphere, espe- 
cially in the centers of galaxy clusters, is still very limited. In preparation for the upcoming 
new generation of radio telescopes (like EVLA, LWA, LOFAR and SKA), a huge effort is 
being made to learn more about cosmological magnetic fields from the observational perspec- 
tive. Here we present the implementation of magneto-hydrodynamics in the cosmological 
SPH code GADGET (Sprin gel et al.ll200ll:ISpringelH2005) . We discuss the details of the im- 
plementation and various schemes to suppress numerical instabilities as well as regularization 
schemes, in the context of cosmological simulations. The performance of the SPH-MHD code 
is demonstrated in various one and two dimensional test problems, which we performed with 
a fully, three dimensional setup to test the code under realistic circumstances. Comparing so- 
lutions obtained using ATHENA (IStone et al J 12008). we find excellent agreement with our 
SPH-MHD implementation. Finally we apply our SPH-MHD implementation to galaxy clus- 
ter formation within a large, cosmological box. Performing a resolution study we demonstrate 
the robustness of the predicted shape of the magnetic field profiles in galaxy clusters, which 
is in good agreement with previous studies. 
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1 INTRODUCTION 

Magnetic fields have been detected in galaxy clusters by radio ob- 
servations, via the Faraday Rotation Signal of the magnetized clus- 
ter atmosphere towards polarized radio sources in or behind clus- 
ters (see Carilli & Taylor 2002; Govoni & Feretti 2004, for recent 
reviews) an d from diffuse synchrotron emission of the cluster atmo- 
sphere (see iGovoni & Ferettill2004l ; iFerrari et alj|2008i . for recent 
reviews). Our understanding of the origin of cosmological magnetic 
fields is particularly limited. But their evolution and possible im- 
plications for structure formation are also not yet fully understood. 
In addition their interplay with the large-scale structure formation 
processes, as well as their link to additional dynamics within the 
cluster atmosphere is unclear, especially their role in the cool core 
regions and the influence of these regions on the evolution of the 
magnetic fields. 

The upcoming, new generation of radio telescopes (like 
EVLA, LWA, LOFAR and SKA) will dramatically increase the 
volume of observational data relevant for our understanding of 
cosmological seed magnetic fields in the near future. To investi- 
gate the general characteristics of the magnetic fields in and be- 
yond galaxy clusters at the level required for a meaningful com- 
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parison to current and forthcoming observations, numerical simu- 
lations are mandatory. Non-radiative simulations of galaxy clus- 
ters within a cosmological environment which follow the evo- 
lution of a primordial magnetic seed field have bee n performed 
using S mooth-Pa rticle-Hydrody namic s (SPH) codes jPolag et al.l 
11999, '2002', ^20 01) as well as Adaptive Mesh Refinern ent (AMR) 
codes (Bruggen 'rtaLl |2005| ; IPubois & TeYssieill2008l) . Although 
these simulations are based on quite different numerical tech- 
niques they show good agreement in the predicted properties of 
the magnetic fields in galaxy clusters. When radiative cooling is 
included, strong amplification of the magnetic fields inside the 
cool-core region of clusters is found (Dubois & TeyssieilboOSl) . in 
good agreement with previous work (Dolag 2000). Cosmological, 
magneto-hydrodynamical simulations were also performed using 
finite-volume and finite-difference methods. Such si mulations are 
used to either follow a primordial magnetic field tLi et al.ll2008l) 
or the creation of magne tic fields in shocks t hrough the so-ca lled 
Biermann battery effect jKulsrud et al]| 19971 : iRvu et al.lll998h . on 
which a subsequent turbulent dynamo may operate. The latter pre- 
dict magneti c field strength in filaments with somewhat higher val- 
ues (e.g. see 'Sigl et al.''2004^ than predicted by simulations which 
start from a primordial magnetic seed field, but are in line with pre - 
dictions of magnetic field values from turbulence IRyu et al.l2008h . 
Therefore further investigations are needed to clarify the structure. 
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evolution and origin of magnetic fields in the largest structures of 
the Universe, their observational signatures, as well as their inter- 
play with other processes acting in galaxy clusters and the large 
scale structure. 

The complexity of galaxy clusters comes principly from their 
hierarchical build up within the large-scale structure of the Uni- 
verse. In order to study their formation it is necessary, to follow a 
large volume of the Universe. However, one must also describe cos- 
mic structures down to relatively small scales, thus spanning 5 to 6 
orders of magnitudes in size. The complexity of the cluster atmo- 
sphere reflects the infall of thousands of smaller objects and their 
subsequent destruction or survival within the cluster potential. Be- 
ing the source of shocks and turbulence, these processes directly 
act on the magnetic field causing re-distribution and amplification. 
Therefore realistic modelling of these processes critically depends 
on the ability of the simulation to resolve and follow correctly this 
dynamics in galaxy clusters. 

Starting from a well-established cosmological n-body 
smoothed particle hydrodynamic (SPH) code GADGET 
dSpringel et alj I2OOII : ISpringell |2003) we present here the im- 
plementation of magneto-hydrodynamics, which allows us to 
explore the full size and dynamical range of state of the art 
cosmological simulations. GADGET also allows us to turn on 
the treatment of many additional physical processes which are of 
interest for structure formation and make interesting links with the 
treatment of magnetic fields for future studies. Th is incl udes ther- 
mal conduction ( Jubelgas et al. 2004; Dolag e t al.ll2004l) . physical 
viscosity ( Sijacki & Springel 2006), cooling and star-formation 
dSpringel & Hernguisli |2003j), detail ed modelling of the stellar 



population and chemical enrichment jTornatore et alj 2004 
and a self cons istent treatment of cosmic rays jEnBlin et alj 



2007) 



2007 



IPfrommer et al.i,2007) . The MHD implementation presented here 
is fully compatible with all these extensions, but here we want 
to focus on non-radiative simulations. All such processes are 
expected to increased the complexity and lead to interplay with 
the evolution of the magnetic field. This would make it impossible 
to critically check the numerical effects caused by the different 
SPH-MHD implementations and therefore we will ignore such 
additional processes in this work. 

The paper is structured as follows: In section 2 we present 
the details of the numerical implementation, whereas in section 
3 we present various code validation tests, all performed in fully 
three dimensional setups. In section 4 we present the formation of 
a galaxy cluster as an example for a cosmological application be- 
fore we present our conclusions in section 5. In addition we present 
a convergence test for the code in the appendix. 



2 SPH-MHD IMPLEMENTATION 

We have implement ed the MHD equations in the cosm ological 
SPH code GADGET jSpringel et al]|200ll : ISpringelll2005l) . In this 
section we present the relevant details of this implementation. 
While developing the MHD implement ation made use of GADGET- 
1 jSpringel et alEoOl.) and GADGET-2 jSpringelll2005) . all simula- 
tions presented in this paper are based on the most recent version 
of the code, GADGET-3 (Springel, in prep). Note that the imple- 
mentation therefore is fully parallelized and benefits from many 
optimizations within the general parts of the code, especially the 
calculation of self gravity and optimization in data structures as 
well as work-load balancing. Therefore, this implementation is an 
ideal tool to follow the evolution of magnetic fields and allows us 



to explore dynamical ranges of more than 5 orders of magnitude 
within cosmological simulations. 

During the last years, many general improvements in the 
implementation of the SPH method have been made. Examples 
include are more modem formulations of the artificial viscosity 
( lMonaghaj|l997h . the introduction of self-consistent correction 
terms from varyi ng smoothing length jSpringel & Hemauistll2002l : 
Monaghan 2002) or the continuous definition of the smoothing 
length ( Spring el & Hernquis t 2002). All these improvements not 
only increased the accuracy and stability of the underlying SPH 
formulation but also improved, directly or indirectly, the accu- 
racy and stability of any MHD implementation. The main, indi- 
rect benefits of these improvements are the self-consistent treat- 
ment of the m agnetic waves within the formulation of the artifi- 
cial viscosity jPrice & Monaghar]|2004al) . the possibility to drop 
the viscosity limiter (see equation 10 and discussion in section 2.3) 
and the inclusion of the correction terms from varying smoothing 
length in both the inductio n equation and the Lorenz force term 
jPrice & Monaghanll2004bh . In the remaining sections we will dis- 
cuss our SPH-MHD implementation in detail. 



2.1 SPH implementation in GADGET 

The basic idea of SPH is to discretize the fluid i n mass ele- 
ments (mi), represented by particles at positions Xi l lLucvlll977l : 
iGingold & Monaghad 1 19771) . To build continuous fluid quanti- 
ties, one starts with a general definition of a kernel smoothing 
metho d. The most frequently used kernel is the B2- 

Spline jMonaghan & Lattanzid[l983) . which can be written as 



W(x,h) 




s; f ^ 0.5 
0.5 < f ^ 1 

1 f 



,(1) 



It is worth stressing that, contrary to other SPH implementation, 
GADGET uses the notation in which the kernel W{x,h) reaches 
zero at a;//i = 1 and not sX x/h — 2. The density pi at each 
particle position Xi can be estimated via 



J 



(2) 



where the smoothing length hi is defined by solving the equation 

(3) 



47r,3 

— hi pi = Nrrii 



A typical value for A'' is in the range of 32-64, which correspond 
to the number of neighbors which are traditionally chosen in SPH 
implementations . 

In GADGET, the equation of motion for the SPH particles 
are implemented based on a derivation from the fluid Lagrangian 
jSpringel & IIernquislll2002l) and take the form 



dvi 
'dt 



(hyd) 



N 



Pt 



The coefficients fi are defined by 



1 + 



hi dpi 
3pi dhi 



(4) 



(5) 



and reflect the full, self-consistent correction terms arising from 
varying the particle smoothing length. The abbreviation Wi — 

W{\fi — fj\, hi) and Wj = VF(|fi — fj\, hj) are the two kernels 
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of the interacting particles. The pressure of each particle is given 
by Pi — AipJ, where the entropic function Ai stays constant for 
each particle in the absence of shocks or other sources of heat. 

To capture shocks properly, artificial viscosity is needed. 
Therefore, in GADGET the viscous force is implemented as 



V At ) 



mjHijViWi- 



(6) 



where Ilij ^ is non-zero only when particles approach each other 
in physical space. The viscosity generates entropy at a rate 



dt 



N 

^' j=i 



(7) 



Here, the symbol Wij denotes the arithmetic mean of the two ker- 
nels Wi and Wj . 

For the parameterization of the artificial viscosity, starting 
with v ersion 2 of GADGET, a formulation proposed by Monaghan 
( Il997h based on an analogy with Riemann solutions of compress- 
ible gas dynamics, is used. In this case, the resulting viscosity term 
can be written as 



n. 



P^3 



- f 

■IK 



(8) 

otherwise, i.e. the pair-wise vis- 
cosity is only non-zero if the particles are approaching each other. 



for fi. 



Here fj.ij 



is the relative velocity projected onto the 



separation vector and the signal velocity is estimated as 



Ci + c, 



(9) 



with Ci = 'jPi I pi denoting the sound velocity. In GADGET-2 the 
values oi = \ and /3 = 3 are commonly used for the dimension- 
less parameters within the artificial viscosity. Here we have also 
included a viscosity-limiter = + ,ff °^')/2, which 

is often used to suppress the viscosity locally in regions of strong 
shear flows, as measured by 



n 



|(V-i;)j + |(Vx^T)j+a. 



(10) 



which can help to avoid spurious angular momentum a nd vortic- 
ity transport in gas disks ( lBalsara|[T99a : ISteinmetz|[T99^ . with the 
common choice Oi = 0.0001ci//ii. 

This also leads naturally to a Courant-like hydrodynamical 
time-step 



(hyd) 



(11) 



where Ccourant is a numerical constant, typically choosen to be in 
the range 0.15 - 0.2. 



2.2 Co-moving variables and integration 

The equations of motion are integrated using a leap-frog integra- 
tion making use of a kick-drift-kick scheme. Within this scheme, 
all the pre-factors due to the cosmological background expansion 
are taken into a ccount within the calculation of the kick- and 
drift-factors (see lOuinn et alj|l997l ; ISpringellbOOSi) . For the inte- 
gration of the entropy within a cosmological simulation, a factor 
[Ha?)^^ ~ si present in equation [7] to take into account that 
the internal time variable in GADGET is the expansion parameter a. 



The formulation of the MHD equations within GADGET has to be 
be adapted accordingly to this choice of variables. 



2.3 Magnetic signal velocity 

A natural generalization of the signal velocity v°^^ in the framework 
of MHD is to replace t he sound velocity a by the fastest magnetic 
wave as suggested by IPrice & Monaghanl ( l2004al) . Therefore the 
sound velocity Ci gets replaced by 



1 

71 



2 I 



POPa 




^ cl{B-n,/\n,\y 



(12) 



As this new definition of the signal velocity also enters the time- 
step criteria i ll It . no extra time-step criteria due to the magnetic 
field has to be defined. We note that we still see improvements in 
the solution to the test problems, if we choose more conservative 
settings within the Courant condition. Therefore we generally use 
Ccourant = 0.075, which is half the value usually used in pure 
hydrodynamical problems. Different authors also propose to use 
different values for a and (3 within the artificial viscosity defini- 
tion l[8j. Whereas typically a = 1 is cho sen, Monaghan 
proposed to use /3 = 3. P rice & Monaghan (2004a.b) propose to 
use [3 — 2 or l3 — 1 respectively. We find slight improvements in 
our test problems when using a — 2 and /3 = 1.5, which we use 
throughout this paper. We also note, that the viscosity suppression 
switch /f/'^'"^ was introduced based on an earlier realization of the 
artificial viscosity and it is not clear if it is still needed. As we note 
significant improvements in our test problems when neglecting this 
switch we do not use this switch throughout this paper. Also for the 
cosmological application presented in the last part of this paper, 
this switch was always turned off. 



2.4 Induction equation 

The evolution of the magnetic field is given by the induction equa- 
tion. 



dB 
lit 



{B ■\/)v~ B{\/ ■ v) 



(13) 



if ohmic dissipation is neglected and the constraint V • B = is 
used. The SPH equivalent reads 



dt 



1 fr 

Ha^ Pi 

N 

^m-j [Bi(vij ■ ViWi) - Vij{Bi ■ ViWi)] 
-2B„ 



(14) 



where (ffa^)~^ = ^ takes into account that the internal time vari- 
able in GADGET is the expansion parameter a. Note that here, by 
construction, only the kernel Wi and its derivative is used. The sec- 
ond term — 2Bi accounts for the dilution of the frozen in magnetic 
field due to cosmic expansion. Both these additions - the {Ha^)~^ 
factor and the —2Bi term - are only present in the cosmological 
simulations and absent for the code evaluation presented in section 
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[3] In component form the induction equation reads 
1 fl 



dBl 
dt 



AT 

mj{vijB\ 



I Vij 



du Wii 



(15) 

Note that, as also suggested by IPrice & Monaghaj ( l2004bh . we 
wrote down the equations including the correction factor ° which 
reflects the correction terms (^^) arising from the variable parti- 
cle smoothing length. Unfortunately it is not possibile to directly 
infer the exact form of the correction factors f rom first principles 
for the induction equation. However. IPrice & Monaghan (2004b) 
showed that, if not chosen in the same way as for the Lorenz force, 
inconsistency between the induction equation and magnetic force 
results. The effect of these factors in the induction equation is quite 
small, but nevertheless one notices tiny improvements in test prob- 
lems when they are included. Therefore, we included them for all 
applications presented in this paper. 



tation of the magnetic force. However their performance was found 
to depend on the details of the simulation setup. In the next sections 
we will briefly discuss the different possibilities in the context of 
building up an implementation for cosmological simulations. 



2.6.1 Adding a constant value 

One method to remove t he instability was pointed out by 
IPhillips & MonaghanI digsj) . who suggested to calculate the maxi- 
mum of the magnetic stress tensor an d to subtract it globa l y from aU 
particles. Or similar, as suggested in Price & Monaghari ( l2005h , to 
subtract the contribution of a constant magnetic field. This is sim- 
ple and straight forward if there is a strong, external magnetic field 
contribution from the initial setup, which can be associated with 
the term one subtracts. However, with cosmological simulations in 
mind, this approach is not very viable, and therefore we did not use 
this approach. 



2.5 Magnetic force 

The magnetic field acts on the gas via the Lorenz force, which can 
be written in a sy mmetric, conservative form involving the mag- 
netic stress tensor jPhillips & MonaghajfTgsl) 



(16) 



The magnetic contribution to the acceleration of the i-th particle 
can therefore be written as 



\ (mag) 



MO 



-2 



\7^W^ 



(17) 



Here needed to transform the equations to the internal 

variables for cosmological simulations and is set to one in all other 
cases. Also /io has to be chosen properly as 

[TIME] 2 [LENGTH] 



47r[MASS]/i2 



(18) 



with h = 1 for non cosmological runs. The factors /™ reflect the 
correction terms (^j^) arising from th e variable particle smoothin g 
length as introduced already (see also lPrice & Monaghan|[2004bh . 
In component form the equation reads 

(vise) 



dt 



Mo 



du 



du \fii\ 



(19) 



It is well known that this formulation becomes unstable 
for situations, in which t he magnetic forces are dominating 
dPhillips & MonaghanlfTgsl) . The reason for this is that the mag- 
netic stress can become negative, leading to the clumping of parti- 
cles. Therefore, some additional measures have to be taken to sup- 
press the onset of this instability. 



2.6.2 Anti clumping term 

iMonaghanI j2000l) suggested the introduction of an additional term 
in the momentum equation which prevents particles from clump- 
ing in the presence of strong magnetic stress. Including this term, 
equation ( 116b reads 



RiBi&i 



(20) 



(21) 



where Ri is a steepened kernel which can be defined as 

The modification of the kernel is made so that contributions are 
significant only at distances below the average particle spacing ui, 
so Wi is defined as Wi = W(ui). Typical values forthe remaining 
parameters are e = 0.8 an d n = 5. This method was also used in 
IPrice & MonagharJ Jiooi^lbl) where they found u\ = 1.5/itobea 
good choice for ID simulations and switched to u\ = l.l h for 2D 
simulations. In agreement with lPrice & MonaghanI ( |2005|) we find 
that in 3D and allowing the smoothing length to vary, this approach 
does not help to suppress the instabilities efficiently. 



2.6.3 div (B) force subtraction 

lB0rve et"ai] ( 1200 ih suggested explicitly subtracting the effect of 
any numerically non-vanishing divergence of B. Therefore, one 
can explicitly subtract the term 



dt 



37 



—PBi m, 
MO 



Pi 



fl 



I Bi 



(22) 



from the momentum equation. Here again, a ~' — and mo are 
introduced to transform the equation to the internal units used. To 
be consistent with the other formulatio ns, we included /™ , which 
are the terms. In the original work te0rve et alj|200lh . /3 = 1 
was choosen. In component form this equation reads 



2.6 Instability corrections 

There are several methods proposed in the literature to suppress the 
onset of the clumping instability which is caused by the implemen- 



dt 



(corr) 



N 



37 



— PBi ^ m, 

Mo 



, Bl dWi rij 



Pt 



du 
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p? 



du 



(23) 



In principle, tliis term breaks thie momentum conserving form of 
the MHD formulation. Howe ver, in practice, this seems to be a mi- 
nor effect. lB0rve et alj j2004h argued that stability for linear waves 
in 2D can be safely reached even when not subtracting the full term 
but choosing /3 < 1; e.g. they suggested j3 — 0.5 to further min- 
imizing the non-conservative contribution. However, it is not clear 
if this stays true for 3D setups and in the non-linear re gime. Addi- 
tional we do not use a higher order kernel as done in lB0rve et al.l 
( |2004|) . therefore it is not clear if their conclusions still hold in our 
case. As the results in our test problems seem unharmed by the pos- 
sible violation of momentum conservation due to the formulation 
of the correction factor in the Lorenz force, we k eep it in the form 
suggest ed in earlier work ( e.g. lB0rveetalJl2OOll) i.e. (3^1. 

InjB^rveetaL I ( l2006l) a more general formalism to obtain /3 
for each particle was introduced with the aim to further minimize 
the violation of momentum conservation in the formulation of the 
correction terms in the Lorenz force. Unfortunately, in the light of 
cosmological simulation, this seems not to be very practical as it 
contains a scan for a maximum value over all particles, which in 
a cosmological context makes no sense as there is not a specific 
single object to which such characteristics can be tuned to. 

In general we find that this correction term significantly im- 
proves all results in our test simulations and effectively suppresses 
the onset of the clumping instability. It was also alre ady success- 
fully used in previous, cosmological applica tions (e.g. lDolag et al.l 
|2004 iRordorf et alj2004lDolag et alj2005h . 



2.7.2 Particle splitting 

lB0rve et"^ ( 1200 1[) developed a scheme to regularize the interac- 
tion of particles in SPH based on a discretization of the smooth- 
ing length hi by factors of 2. In such cases, interactions between 
two particles with different smoothing length can be realized by 
splitting the one with the larger smoothing length into 2" (where 
n is the dimensionality) particles, placed on a 2" sub-grid. Such 
split particles then have the same smoothing length as the particles 
with which they interact. Originally this scheme was invented to 
avoid the problems induced by a variable smoothing length (before 
the correction terms where prope rly introduced in SPH) and gav e 
good results in ID and 2D (see lB0rve et"Zll200ll. |2004 l2006h . 
However, in 3D the resulting change in the number of neighbors 
is quite large when the smoothing length is quantized in factors 
of 2. Therefore the additional sampling noise for particles can be 
large. This is especially problematic when the lower density is ap- 
proaching - but still above - the threshold for doubling the particle 
smoothing lengths. Here particles have a particularly small smooth- 
ing length (relative to the optimal, unquantized one) and therefore 
have only small numbers of neighbors. Unfortunately, this effect is 
much larger than the gain in accuracy by the r egularization, at least 
when based on standard SPH formalism, (see lDel Pral2003l) . Also, 
as the terms formally take care of all correction terms induced 
in the formalism when allowing a variable smoothing length, this 
splitting - and specifically the quantization of the smoothing length 
- is no longer needed. This might be different when further improv- 
ing the SPH (and specially the MHD) method. For example when 
re-mapping techn iqes based on Voronoi tessellation are used (e.g. 
lB0rveetalJl2OO6^ or special c oordinates like spherical or cylindi- 
cal are used (e.g. Omang et al.ll2006) . 



2.7 Regularization schemes 

Beside instabilities, noise (e.g. fluctuations of the magnetic field 
imprinted by numerical effects when integrating the induction 
equation) is a source of errors in SPH-MHD implementations. The 
goal of a regularization scheme is to obtain a magnetic field which 
does not show strong fluctuations below the smoothing length. This 
can be achieved indirectly through improvements in the underlying 
SPH formalism and by reformulation of the the interactions to re- 
duce the creation of small irregularities from numerical effects (e.g. 
through particle splitting). Alternative approaches are to directly 
suppress the magnetic field or to dissipate small irregularities by 
introducing artificial dissipation. 



2. 7. 1 Improvements in the underlying SPH formalism 

Here , the entropy conserving formalism (Sprin gel & Hemguisll 
l2002h of the underlying SPH implementation contributes to a sig- 
nificant improvement of the MHD formalism compared to previous 
MHD implementations in SPH by generally improving the den- 
sity estimate and the calculation of derivatives. It has to be noted, 
that this is not only due to the terms, but in large part also 
from the new formalism for calculation of the smoothing length. 
As described before, the smoothing length hi for each particle is 
no longer calculated by counting neighbors within the sphere, but 
by solving equation Jsj for each formal number of neighbors N, 
there exists only one unambiguous value of hi. Note that, as this 
equation is solved iteratively, it is usual to give some allowed range 
of A'^, however in our case we can choose the range smaller than 1 
and typically we use = 64 ± 0.1. 



2.7.3 Smoothing the magnetic field 

Another method to remove small scale fluctuations and to regular- 
ize the magneti c field is to smooth the magnetic field periodicly. 
As suggested bvl B0rve et al.l ( 1200 ll) . one can calculate a smoothed 
magnetic field (Si) for each particle. 



(24) 



Then, in periodic intervals, one can calculate a new, regularized 
magnetic field by 



Br" = <?(B«) + (i-g)B.. 



(25) 



Note that this, in principal, acts similar to the mixing process on 
resolution scale present in Eulerian schemes. However, introduced 
in this way, the amount of mixing (e.g. dissipation) of magnetic 
field depends on the frequency with which this procedure is ap- 
plied and the value of q chosen. Typically, we set q to one and per- 
form the smoothing at every 15*''-20*'' main time-step. It is worth 
while to mention that implemented in this form, total energy is not 
conserved (as magnetic field fluctuations on scales smaller than the 
smoothing length are just removed) and, as the time-steps depend 
on the chosen resolution, this method is even resolution dependent. 
Never the less it leads to improvements in the results of our test 
problems, without strongly smoothing sharp features. It also works 
without pro blem in 3D and has alrea dy been used in cosmological 
simulations dDolag et al.ll2004l2005h . 
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Figure 1. T est 5A at time t = 7 with the MHD implementation similar to 
that used in lDolagetaljfT999ll200A . but already including the instability 
correction due to subtraction of the div(_B) term in the force equation and in 
a fully three-dimensional setup. Shown in the first row are the density (left 
panel), total energy and pressure (middle panel) and the x,z component of 
the velocity field (right panel). The second row shows the i/-component in 
the velocity field (left panel), the three components of the magnetic field 
(middle panel) and the measure of the div(B) error, see equation i32\ . in 
the right panel. The black lines with error bars show the SPH results, the 
red lines are the reference results obtained with Athena in a ID setup. 



2.7.4 Artificial magnetic dissipation 

A nother possibility to regulari ze the magnetic field was presented 
bv lPrice & MonaghanI ( l2004al) . who suggested including an artifi- 
cial dissipation for the magnetic field, analogous to the artificial vis- 
cosity used in SPH. In lPrice & MonaghanI ( l2004a[) it was suggested 
that the dissipation terms be constructed based on the magnetic 
field component perpendicular to the line joining the interacting 
particles. However, to better suppress the small scale fluctuations 
within the magnetic field which app ear due to numerica l effects 
especially in multi-dimensional tests, IPrice & MonaghanI ( l2004bh 
suggested basing the artificial dissipation on the change of the total 
magnetic field rather than on the perpendicular field components 
only. We also found this to work significantly better in our test 
cases and therefore only use the later implementation throughout 
this paper. Such an artificial dissipation term can be included in the 
induction equation as 



At 



(diss) 



1 PiOlB 



{B,-Bj)-^- V.m. (26) 



The parameter as is used to control the strength of the effect, typ- 
ical values are suggested to be around q_b ~ 0.5. Similar to the 
artificial viscosity, this will create entropy at the rate 

(diss) 
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\ At 
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Pi, 
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The pre-factor (7 — 1)/(p7~^) properly converts the dissipation 
term to a change in entropy. 

This method reduces noise significantly. However, depending 
on the choice of as, it can also lead to smearing of sharp features. 



Figure 2. As figure [T] but including the magnetic waves in the signal ve- 
locity and turning off the shear viscosity suppression as explained in the 
code description. The main advantage is a significant reduction in the noise, 
specifically in the velocity, but also in the magnetic field. Also the div(i?) 
errors are reduced by a factor of 2. 



To avoid this outside of st rong shocks (e.g. where this is needed), 
IPrice & MonaghanI ( l2005h proposed evolving as for each parti- 
cle, simil ar to the hand ling of the time dependent viscosity as sug- 
gested by iMorris & M onaghan ( 1997). Such, evolution of as for 
each particle will be followed by integrating 



das 



(as 



+ S, 



where the source term 5* can be chosen as 



S = So max 



(28) 



(29) 



(see IPrice & Monaghanll2005l) . The time-scale r defines how fast 
the dissipation constant decays. Taking the signal velocity, one can 
translate this directly into a distance to the shock over which the 
dissipation constant decays. A useful choice of r can be written as 

hi 



C Usi 



(30) 



where the constant C typically is chosen to be around 0.2, allowing 
the dissipation constant to decay within a ti me that corresponds 
to the shock travelling 5 kernel lengths (see IPrice & MonaghanI 
l2004ah . 



2.8 Euler potential 

A very elegant way to implement the MHD equations in La- 
grangian codes is the usage of so called Euler potentials (see 
.Rosswog & Price 2007. , and references therein). Two independent 
variables a and /3 are constructed to correspond to an implicit 
gauge for the vector potential. They can be thought of as labels 
of magnetic field lines and will be adverted with the flow. In this 
formulation, the magnetic field at any time can be represented as 



B = Va X V/3. 



(31) 



In principle, having obtained the magnetic field, one could also use 
this magnetic field in the equation of motion as before. However, 
this would mean that the magnetic force is based on the second 
derivative of a variable. This is usually quite noisy and not recom- 
menda ble unless regularization schemes are implemented as done 
by e.g. IRosswos & Pried ( l2007h . In addition, due to its form, an 
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Figure 3. As figure |2] but including the regular smoothing of the magnetic 
field as a regularization scheme. Th i s SPH-MHD implementation basically 
reflects the one used in iDolag et alj j2004l200a . The main advantages are 
a further, significant reduction in the noise as well as a strong reduction of 
the div(_B) eiTors by a factor of 10. 



Figure 4. As figure|2] but including artificial magnetic dissipation as a reg- 
ularization scheme. Similar to the smoothing of the magnetic field, signif- 
icant reduction in the noise as well as a strong reduction of the div(iJ) 
errors by a factor of 10 is obtained compared to the basic SPH-MHD 
implementation. 



implementation based on the Euler potential cannot be easily inves- 
tigated in 3D test problems with periodic boundary conditions for 
the resulting magnetic field as can be done for the other implemen- 
tations. A simple example here is constant magnetic field, which 
can be represented by a linear function of the Euler potentials. 

Therefore we use this simple description only as a check in 
cosmological simulations, to investigate the influence of div(i3) 
driven errors. Applied to cosmological simulations the evolution of 
the magnetic field predicted when using Euler potentialsis an upper 
bound on the amplification processes in the absence of any dynamo 
action. Therefore Euler potentials are a useful tool to check the 
influence of numerics on the results of cosmological simulations 
where we have no other means to verify the results. 



3 TEST PROBLEMS 

To test performance of the code and to infer the optimal nu- 
merical settings for the regularization schemes, we performed 
the se ries of shock-tube problems as presented by IRvu & Jones 
1995t) . In particular test JA, which is also used in iBrio & Wi] 
1988h . was used to show the effects of different numerical treat- 



ments. Additionally we performed several 2D test cases includ- 
ing the F ast Rotor test (Toth 2000, ; iLondrillo & Del Zannj|2OO0l: 
Balsara & Spicerl Il999i), a Stro ng Blast jLondrillo & Del Zannal 
'200a IBa lsara & Spicer' 'l999^ and the Orzang-Jang_Vortex 
jOrzang & Tang 1979; Dai & Woodward 1994; Pico ne & Dahlburd 



I1991I ; iLondnllo & Del Za nna 2000). To obtain results under realis- 
tic circumstances, we performed all the tests by setting up a fully 
three-dimensional particle distrib ution. We also avoid starting from 
regular grids but used glass-like j White! 0*9961) initial particle dis- 
tributions instead. To obtain such a configuration, the particles are 
originally distributed in an random fashion within the volume and 
then allowed to relax until they settle in a equilibrium distribution 
which is quasi force-free and homogeneous in density. This is simi- 
lar to the distribution of atoms in an amorphous structure like glass. 
Compared to a distribution of the particles based on a grid, this 
guarantees that all the kernel averages in the SPH formalism sam- 
ple the kernel in a uniform way rather than multiple times at the 
same distances (which furthermore would be fractions of the under- 
lying grid spacing). For all tests we used the same particle masses. 



independent of the initial density. Therefore, typical initial particle 
distributions for the shock-tube tests where based on 5^ particles 
in low density and 10"' particles in high density regions within unit 
volume. Usually, these unit volumes are then replicated 35 times 
along the a;-direction each. For some test cases with strong (and 
therefore fast) shocks, we evolved the simulations longer In such 
cases we doubled the simulation setup size in the x-direction. 

We assume ideal gas (e.g. 7 = 5/3) and, as described before, 
use an equivalent of 64 neighbors for calculating the SPH smooth- 
ing length. This ensures that, in the low density regions, the SPH 
particles get smoothed over a region corresponding to a unit length. 
The number of resolution elements corresponding to a unit length 
therefore ranges from 1 to 4, depending whether one associates the 
smoothed region or the mean inter-particle distance with the effec- 
tive resolution in SPH. In general, SPH converges somewhat slower 
compared to grid codes when comparing simulations with the same 
number of grid cells as SPH particles (see Appendix A for an ex- 
ample). 

For the SPH results we usually plot the mean within a 3D slab 
corresponding to the smoothing length and (as error bars) the RMS 
over the individual particles withi n this volume. Th e reference so- 
lution was obtained using Athena jStone et^l2008l) with typically 
10-20 resolution elements per unit length, depending on the indi- 
vidual test. As one criteria of the goodness of the SPH simulation 
result we use the usual measure for the non-vanishing divergence 
of the magnetic field, 

i5^^ = div(B)A. (32) 



3.1 Shock tube 5A 

The most common ly used MHD sh ock-tube test is the o ne used by 
iBrio & Wj l ll988h . e.g. test 5A in IRvu & JonesI ( Il995h . The rea- 
son for this is that it involves a shock and a rarefraction moving 
together. Therefore it allows simultaneous testing of the code in 
different regimes. 

Figure [T] shows the result for a code implementation similar 
to the first implementation used to study galaxy clusters (e.g. see 
iDolag et"ar . 1999, 2002). In addition, the instability correction due 
to subtraction of the div(i3) term was used in the force equation. 



8 K. Dolag, F. Stasyszyn 




-20 -10 *0 20 30 -20 -"0 tO 20 30 -20 -10 10 20 30 



Figure 5. As figure|2] but including time dependent artificial magnetic dissi- 
pation as a regularization scheme. No significance improvement is obtained. 
Note that here in the lower right panel the artificial dissipation constant 
{as ) is shown. The effect of suppressing the dissipation is clearly visible, 
and the maximum value is only reached in peaks associated with the re- 
gion of strong shocks. However the improvement in the smearing of shaip 
features is not very significant. 



Various hydro-dynamical variables at the final time (e.g. t — 7 
in this case) are shown. The black lines with error bars show the 
SPH-MHD result, the red lines are the reference result obtained 
with Athena in a ID setup. Shown are (from upper left to the bot- 
tom right panel) the density, total energy and pressure, the x- and 
z-component of the velocity field, the j/-component in the veloc- 
ity field, the three components of the magnetic field and the mea- 
sure of the div(i3) error, obtained from equation ([32}. Here we 
also switched back to the conv entional formulati on of the artificial 
viscosity as described by e.g. iMonaghanI jl992h . rather than that 
based on signal velocity as used in GADGET-2. Although the SPH- 
MHD results in general follow the solution obtained with Athena, 
there is a large scatter in the individual particle values within the 
3D volume elements, as well as some instability, especially in the 
low-density part. But note that although the mean values for the 
internal energy, as well as the velocity or magnetic field, can lo- 
cally show some systematic deviations from the ideal solution, the 
total energy shows much better, nearly unbiased, behaviour. This 
demonstrates the conservative nature of the symmetric formula- 
tions in SPH-MHD. 

Noticeable reduction of noise is obtained when using the 
signal-velocity based artificial viscosity and including the magnetic 
waves in the calculation of the signal velocity. Therefore, the mag- 
netic waves are directly captured for the time step calculation and 
in the artificial viscosity, needed to capture shocks. Also switch- 
ing off the shear viscosity suppression again leads to significant 
reduction in scatter. This can be seen in figure [2] where the noise 
in the velocity as well as in the magnetic field components is sig- 
nificantly reduced. Values of div(i3) are also reduced (by a factor 
of « 2) compared to before. In general, the SPH-MHD implemen- 
tation gains from the new formulation of SPH, including the 
terms and the new way to determine the SPH smoothing length, 
both contributing to a reduction of noise (and div(i?)) in the gen- 
eral treatment of hydro-dynamics. We will refer to this implemen- 
tation of SPH-MHD as basic SPH-MHD further on in this paper. 



3.2 The effect of regularization 

As described in section l l2.7t . there are several suggestions for reg- 
ularization of the magnetic field. Here we will show results ob- 
tained by two regularization methods, namely smoothing the mag- 
netic field in regular intervals and including an artificial dissipation. 

For the first method, the magnetic field is smoothed using 
the same kernel as used for the normal SPH calculations. In this 
case, there are two numerical parameters one can choose. One is 
q in equation l |25t , which quantifies the weight with which the 
smoothed component enters into the updated magnetic field. We al- 
ways use g = 1 here, which means that we completely replace the 
magnetic field by the smoothed value. The second is Tbs, which 
is the time interval at which the smoothing is done. Here we use 
a value corresponding to a smoothing every 30*'' global time step. 
This correspond to the SPH-MHD implementation used to study 
the magnetic field in clusters and large scal e structure within the 
local universe, see iDolag etal.ll2004l2005l) . Figure [3] shows the 
result for the same shock-tube test as before. Clearly, the noise 
in the individual quantities is strongly reduced. Also the error in 
div(i3) is reduced by more than one order of magnitude. Note that 
the error bars for the SPH-MHD implementation are of the size 
of the line width or smaller in most of the cases and therefore no 
longer clearly visible. However, one can notice some small effect 
of smearing sharp features. Additinalally, some states - like the re- 
gion with the negative x-component of the velocity behind the the 
fast rarefaction wave propagating to the right - converge to values 
which have small but systematic deviations from the exact solution. 

In the second method, the magnetic field can be dissipated in 
the same way as artificial dissipation works in the hydrodynamics. 
Here the numerical parameter one has to chose is the strength of this 
artificial, magnetic dissipation as in equation i26\ and l |27| l. Fig- 
ure|4]shows the result for the same shock-tube test as before using 
ub ~ 0.1. Similar to the first regularization method presented, the 
noise in the individual quantities is strongly reduced and also the 
error in div(_B) is reduced by one order of magnitude. Again, the 
error bars are smaller than the line width nearly everywhere. Also, 
some small effects of smearing sharp features are visible as well 
as some small but systematic deviations from the exact solution. In 
general, this method works slightly better than the smoothing of the 
magnetic field, but the differences are generally small. 

One idea to reduce the unwanted side effects of such regular- 
ization schemes was presented in , Price & Monaghan (.2005) and is 
based on a modification of the artificial, magnetic dissipation con- 
stant as. Whereby every particle evolves its own numerical con- 
stant, so that this value can decay where it is not needed and there- 
fore the effects of the artificial dissipation are suppressed here. Fig- 
ure|5]shows the same test as before, but this time where as evolves 
for each particle, as shown in the lower right panel. Clearly, the val- 
ues are strongly reduced outside the regions associated with sharp 
features (e.g. shocks), but the effect of smearing sharp features and 
the small offset of some states are not significantly reduced. This 
is because in the region in which these side effects originate, the 
dissipation is still working with its maximum numerical value. On 
the other hand, due to the suppression of the artificial magnetic dis- 
sipation constant outside the shock region, the regularization after 
the shock passes is nearly switched off. Therefore it is not as effi- 
cient as before in the post shock region, visible as increase in the 
div(i3) error compare to the run with a constant, artificial magnetic 
dissipation. 
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TEST Nr. 






Left 








Right 




P 


V 


o 


r 


P 


V 


ti 




— lA — 


1.00 


[10.0,0.0,0.0] 


[5.0, 5.0, 0.0]/(47r)2 


20.0 


1.000 


[-10.0,0.0,0.0] 


[5.0, 5.0, 0.0]/{47r)2 


1.00 


— lB — 


1.00 


[0.0,0.0,0.0] 


[3.0,5.0,0.0]/(47r)2 


1.0 


0.100 


[0.0,0.0,0.0] 


[3.0,2.0,0.0]/(47r)2 


10.0 


— 2A— 


1.08 


[1.2,0.01,0.5] 


[2.0,3.6,2.0]/(47r)2 


0.95 


1.000 


[0.0,0.0,0.0] 


[2.0,4.0, 2. 0]/{47r)2 


1.00 


— 2B — 


1.00 


[0.0,0.0,0.0] 


[3.0,6.0,0.0]/(47r)2 


1.0 


0.100 


[0.0,2.0,1.0] 


[3.0, 1.0,0.0]/{47r)2 


10.0 


— 3A— 


1.00 


[50.0,0.0,0.0] 


-[0.0, 1.0, 2.0]/(47r)2 


0.4 


0.100 


[0.0,0.0,0.0] 


[0.0, 1.0, 2.0]/{47r)2 


0.20 


— 3B — 


0.10 


[-1.0,0.0,0.0] 


[0.0,1.0,0.0] 


1.0 


1.000 


[1.0,0.0,0.0] 


[0.0,1.0,0.0] 


1.00 


^4A— 


1.00 


[0.00,0.0,0.0] 


[1.0,1.0,0.0] 


1.0 


0.200 


[0.0,0.0,0.0] 


[1.0,0.0,0.0] 


0.10 


— 4B— 


0.40 


[-0.669,0.986,0.0] 


[1.3,0.0025293,0.0] 


0.5247 


1.000 


[0.0,0.0,0.0] 


[1.3,1.0,0.0] 


1.00 


— 4C — 


0.65 


[0.667, -0.257, 0.0] 


[0.75,0.55,0.0] 


0.50 


1.000 


[0.4,-0.94,0.0] 


[0.75,0.00001,0.0] 


0.75 


— 4D — 


1.00 


[0.0,0.0,0.0] 


[7.0,0.001,0.0] 


1.0 


0.300 


[0.0,0.0,0.0] 


[7.0,1.0,0.0] 


0.20 


— Brio Wu — 


1.00 


[0.0,0.0,0.0] 


[0.75, 1.0,0.0] 


1.0 


0.125 


[0.0,0.0,0.0] 


[0.75,-1.0,0.0] 


0.10 



Table 1. Summary table with the initial conditions of the left and right side of the shoclc tubes. 



3.3 Shock tube problems 

As can be seen in figures[3]and|4l the side effects of smoothing fea- 
tures by the different regularization methods depend on the details 
of the underlying structure of the shock-tube test. Even more inter- 
esting, the states where one can see small deviations from the ideal 
solution are different for the two different regularization methods. 
Therefore w e performed the full set of different shock-tube tests as 
presented in lRvu & Jonej jl995 ^ to test the overall performance of 
the different implementations under different circumstances. The 
four test families deal with different complexities of velocity and 
magnetic field structures, leading to different kinds of waves prop- 
agating. A summary of the results of these tests can be found in fig- 
ure |6] Plotted are the total energy (left panels), the velocity along 
the x-direction (middle panels) and the magnetic field along the 
y-direction (right panels). The red lines reflects the ideal solution 
obtained with Athena, the black lines with error bars mark the re- 
sults from the SPH-MHD implementation using the magnetic field 
smoothing every 30th main time step. Note that the error bars in 
most cases are smaller than the line width. The initial setups for the 
shock-tube tests can be found in table[T] which lists the state vector 
of the left and right states for the different shock tube tests. 

The first family of tests (lA/lB) has no structure in the tan- 
gential direction of the propagating shocks in magnetic field and 
velocity, e.g. = Vz = 0. As we expect, in the lA test, the strong 
shock (large jump in v^) leads to some visible noise in the magnetic 
field component By, also translating into significant noise in the to- 
tal energy. The regularization method here suppresses the formation 
of the intermediate state in By in the SPH-MHD implementation, 
as can be seen in figure [6(a)l The second case, the IB test, the weak 
shock is captured well. Again in some regions some smearing of 
sharp features due to the regularization method is clearly visible. 

The second class of shocks (2A/2B) involve three dimensional 
velocity structures, where the plane of the magnetic field rotates. 
All features (e.g. fast/slow shocks, rot ational discontinuity and 
fast/slow rarefaction wave, for details see lRvu& Jones! ( fl995l) ). are 
well captured, see figure [6(c)l and [6(d)l Some of the features are 
clearly smoothed by the regularization method. 

The third class of tests (3A/3B) shows handling of magne- 
tosonic structures. The first has a pair of magnetosonic shocks with 
zero parallel field and the second are magnetosonic rarefractions. 
Although there is slightly more noise present, all states are captured 
extremely well, except the numerical feature left at the position di- 
viding the two states initially, see figure [6(e)] and |6(f)l 

The fourth test family (4A/4B/4C/4D) deals with the so-called 



switch-on and switch-off structures. The tangential magnetic field 
turns on in the region behind switch-on fast shocks and switch-on 
slow rarefractions. Conversely, in the switch-off slow shocks and 
switch-off fast rarefractions the tangential magnetic field turns off. 
Again, all structures are captured well with the exception of one 
feature in figure [6(h)] as well and maybe [6Q)] too, where clearly the 
regularization leads to the washing out of a state. Otherwise the 
regularization leads to smoothing of some structures similar to the 
tests presented before. 

In general, figure[6]demonstrates that all these different situa- 
tions have to be included when trying to measure the performance 
and quality of different implementations of regularization methods. 



3.4 Finding optimal numerical parameters 

To optimize, we performed all these 1 1 shock-tube tests with vari- 
ous different settings for the parameters in the regularization meth- 
ods and evaluated the quality of the result obtained with the SPH- 
MHD implementation. To measure this, we used two estimators. 
First, we have chosen the mean of all div(i3) errors within the sim- 
ulation region shown in the plots, as defined by 



^div(S) 



B) = (div(B)A 



B\ 



(33) 



Second, we measured the discrepancy of the SPH-MHD result for 
the magnetic field relative to the results obtain by Athena. There- 
fore we calculate first 



B\ 



aa(a;))' 



RMS^, {x) 



(34) 



for each component i of the magnetic field B within each 3D slab 
corresponding to the smoothing length. The RMS of reflects the 
noise of B, within the chosen slab. We then calculate 



3- E^. 




^RMS 



(35) 



for each component of the magnetic field. This includes both con- 
tributions, the deviation of the SPH-MHD from the ideal solution as 
well as the noise within each 3D slab of the SPH-MHD implemen- 
tation. To judge the improvement of the regulatization methods we 
sum up all three components and further relate this measurement to 
the value obtained with the basic SPH-MHD implementation, e.g. 
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Figure 6. Representative plots of the additional 10 shock-tube tests from lRvu & Jones! il99 5l). Shown for each test are the total energy (left panels), the velocity 
along the x-direction (middle panels) and the magnetic field along the j/-direction (right panels). 



As = . - 1. 

We will use these two error estimators, A 
sure the quality of the individual SPH-MHD implementations. 



(36) 



div(fl) 3"'! ^B, to mea- 



3.4.1 Regularization by smoothing the magnetic field 

Choosing the time interval between smoothing the magnetic field 
is a compromise between reducing the noise in the magnetic field 
components, as well as div(B)) (by smoothing more often) and 
preventing sharp features from being smeared out. The left two 
panels of Figure [7] show a summary of the results of the individ- 
ual shock- tube test computed with different smoothing intervals. 
As expected, when using shorter smoothing intervals the error in 



div(_B) reduces. For the quality measure of the SPH-MHD imple- 
mentation the situation changes. Short smoothing intervals gener- 
ally increase the discrepancy, many of them even to larger values 
than the basic SPH-MHD run. Specifically 4B and 4C show strong 
deviations due to smearing of sharp features. Note that the non 
monotonic behavior shown in some tests usually relates to some 
residual resonances between the magnetic waves and the smooth- 
ing intervals in the noise. Some tests show a minimum in the differ- 
ences at smoothing intervals around 20. The test 3A seems to prefer 
even shorter smoothing intervals. In general is not clear how an op- 
timal decision between such quality measure and the reduction in 
div(_B) can be reached, given the different nature and amplitude 
of the two measures. However, ignoring 4B which strongly suffers 
from smearing sharp features when smoothing the magnetic field, a 
good compromise seems to be for values around 20-30, where pro 
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and con in the quality measure are small and compensating within 
the different tests but div(_B) is still drastically reduced in all tests. 
We will refer to this as the Bsmooth SPH-MHD implementation in 
the rest of the paper. 



3.4.2 Artificial dissipation 

As before, choosing the value for the artificial magnetic dissipa- 
tion constant as is a compromise between reducing the noise in 
the magnetic field components (as well as reducing div(_B)) and 
preventing sharp features from smearing out due to the effect of the 
dissipation. The right two panels of Figure |7] show a summary of 
the results of the individual shock-tube tests computed with differ- 
ent values for the artificial magnetic dissipation. As expected, us- 
ing larger values reduces the error in div(i3) significantly. Similar 
to before, using larger values also generally results in an increase 
of the discrepancy between the SPH MHD implementation and the 
true solution, again usually to even larger values than in the basic 
SPH-MHD run. As before, especially the shock-tube test 4B and 
4C show strong deviations due to smearing of sharp features. Note 
that here less non-monotonic behavior is visible (except for test 
4B). The main reason is that dissipation is a continuous process, so 
resonances between dissipation and the magnetic waves cannot be 
very pronounced. As before, it is difficult to infer the best choice 
of parameter for all the tests. Again, once ignoring 4B, a compro- 
mise for choosing ub seems to be between 0.02 and 0.1. Choosing 
a B close to the upper value of 0. 1 might lead to a significant re- 
duction in div(i3) without to strong signature from smearing out 
sharp features. We will refer to this as the dissipation SPH-MHD 
implementation in the rest of the paper. 



3.4.3 Time dependent artificial dissipation 

One idea to reduce the effect of the artificial dissipation is to make 
the artificial magnetic dissipation constant aB time dependent. The 
idea here is that, if the evolution of as is properly controlled, dissi- 
pation will happen only at the places where it is needed and it will 
be suppressed in all other parts of the simulation volume. The evo- 
lution of as is controlled by the two parameters So (source term) 
and C (decay term) where we have chosen as"^ and as''" as 0.01 
and 0.5 respectively. Figure[8]shows the result for varying these two 
parameters. As before, generally, the larger the dissipation is (e.g. 
large source term or small decay time) the smaller the noise and 
the error in div(i?) becomes. However, as soon as these parameters 
have values which drive as in the shocks to the maximum allowed 
value, there is marginally no gain in quality, although the values for 
a s outside the shocks can still be quite small. Therefore, the time 
dependent method does not improve the results significantly, as the 
regions in which the artificial dissipation constant is suppressed do 
not significantly contribute to the smearing of sharp features. 



3.5 Multi dimensional Tests - Planar Tests 

Besides the one dimensional shock tube test described in the pre- 
vious section, two dimensional (e.g. planar) test problems are a 
good test-bed to check code performance. Such higher dimensional 
tests include additional interaction between the evolving compo- 
nents with non-trivial solution. These can be quite complex (with 
several classes of waves propagating in several directions) such as 
the Orszang-Tang Vortex or simple (but with strong MHD discon- 
tinuities) such as Strong Blast or Fast Rotor. 
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Figure 10. Diagonal {x = y) cut through the Fast Rotor at t = 0.1 show- 
ing the density. In black result obtained with ATHENA. The pink line with 
the red error bars show the GADGET solution using the basic SPH-MHD 
implementation. Also here the red error bars reflect the dispersion of the 
values among individual particles within a slab corresponding to the lo- 
cal smoothing length. In general, the SPH-MHD result shows an excellent 
agreement in all the features (peaks, valleys and edges), however there is a 
visible over- smoothing at the outer edges in the GADGET result. 



.J.J. 7 Fast Rotor 



This test problem was introduced by iBalsara & Spiceil ( Il999l) . to 
study star formation scenarios, in particular the strong torsional 
Alfven waves, and is also commonly used to validate M HD im- 
plem entations (for example see Toth 2000; Londrillo & De l Zann j 
Um^ IPrice & MonaghanI l2005l : lB0rve et al.ll200 6). The test con- 
sists of a fast rotating dense disk embedded in a low density, static 
and uniform media, with a initial constant magnetic field along the 
x-direction (e.g. Bx ~ 2.^%~^^'^). In the initial conditions, the disk 
with radius r = 0.1, density p = 10 and pressure P = 1 is spin- 
ning with an angular velocity u — 20. It is embedded in a uniform 
background with p = P = 1. Again we setup the initial conditions 
by distributing the particles on a glass like distribution in 3D, using 
700 X 700 X 5 particles and periodic boundaries in all directions for 
the background particles. The disk is created by removing all parti- 
cles which fall inside the radius of the disk and replacing this space 
with a denser representation of particles of the same mass. As an 
ideal solution to compare with, we again used the result of a sim- 
ple, two dimensional ATHENA run with 400 x 400 cells. A visual 
impression of the results can be obtained from the maps presented 
in figure |9] 

Figure [lO]presents another quantitative comparison. Shown is 
a diagonal cut through the Fast Rotor at f = 0.1 showing the den- 
sity. The different lines show the result obtained with ATHENA 
(black line) and for the three different SPH-MHD implementations 
in GADGET (colored lines). The very small, red error bars reflect 
the RMS of the values held by the individual particles within the 
3D slab through the three dimensional simulations corresponding 
to the local smoothing length. The results show remarkable agree- 
ment between the simulati ons and also compare well w ith results 
quoted in the literature (e.g. lLondrillo & Del Zannal2000l) . Here the 
smoothing of sharp features in the two implementations with regu- 
larization is quite clear visible and leads to a less good match than 
the basic SPH-MHD implementation. 

Note that although we perform our calculations in three di- 
mensions and without a regularization scheme, the implementation 
produce a result, which has the same quality as other schemes in 
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Figure 7. Shown in the two left panels are the mean en'or in divergence (first panel on the left) and the measure of the quality (second panel from left) as 
defined in equation j36) obtained by the SPH-MHD implementation for different values of the smoothing interval. The different lines are for the 1 1 different 
shock-tube tests as indicated by the labels. The two right panels show the same quantities but for different values of the artificial magnetic dissipation constant 
Q-B- 




Figure 8. Similar to figure[2]but for different values of the source term So (left panels) and the decay term C (right panel) of the time dependent, artificial 
magnetic dissipation. 



two dimensions w ith regularization (e.g. IPrice & Monaghanll2005l : 
lB0rve et alj200^ . 



3.5.2 Strong Blast 

The Strong Blast test consists of the explosion of a circular hot 
gas in a static magnetized medium and is also regularly used 
for MHD code validation (see for example iLondrillo & Del Zannal 
Balsara & Spicer 1999). The initial conditions consist of a 
constant density p = 1 where a hot disk of radius ro = 0.125 is 
embedded, which is a hundred times over-pressured, e.g the pres- 
sure in the disk is set to Pd = 100 whereas the pressure outside 
the disk is set to Po = 1. In addition there is initially an overall 
homogeneous magnetic field in the a;-direction, with a strength of 
Bx ~ 10. The system is evolved until time t = 0.02 and an outgo- 
ing shock wave is visible which, due to the presence of the magnetic 
field, is no longer spherical but propagates preferentially along the 
field lines. Figure [TT] shows the density at the final time, compar- 
ing the ATHENA results with the results from the three different 
SPH-MHD implementations in GADGET. Although the setup is a 
strong blast wave, there is no visible difference of the SPH-MHD 
implementation with the ATHENA results. This is quantitatively 
confirmed in figure[T2] which shows a horizontal cut (at y — 0.5) of 



the density through the Strong Blast test, comparing the ATHENA 
(black line) with the GADGET (colored lines with error bars) re- 
sults. Besides very small variations there is no significant difference 
between the two results and all features are well reproduced by the 
SPH-MHD implementation. Note that the error bars of the GAD- 
GET results again are almost in all cases smaller than the shown 
line width. 



3.5.3 Orszang-Tang Vortex 

This planar test problem, introduced by lOrzang & Tan3 ( Il979t) . 
is well known to study the interaction between several classes 
of shock waves (at different velocities) and the transition to 
MHD turbulence. Also, this test is commonly used to vali- 
date MHD implementations (for example see Dai & Woodward 
1994 IPicone & Dahlbiir^ Il99ll: ILondrillo & Del Zanna 2000l : 
Price & MonaghaiJ l2OO5l : (B0rve et alj|2006l) . The initial conditions 
for an ideal gas with 7 = 5/3 are constructed within a unit-length 
domain (e.g. x = [0, l],y = [0, 1]) with periodic boundary con- 
ditions. The velocity field is defined hy — — sm{2ny) and 
Vy = sin(27ra;). The initial magnetic field is set to be Bx = BqVx 
and By 



Bosin^A-Kx) 
the pressure is set to P 



The initial density is p = 7P and 
- jBq. This system is evolved until 
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Figure 9. The magnetic pressure (_B^/2) in the Fast Rotor test at t = 0.1. The ATHENA solution of the test problem is shown in the upper left panel. The 
upper right panel shows the results obtained with the basic SPH-MHD implementation. The lower left and right panels shows the result obtained with the 
Bsmooth SPH-MHD and the dissipation SPH-MHD implementation respectively. All the main features are well reproduced in the GADGET runs. The shape, 
positions and amplitudes correspond quite well, although the GADGET runs appear slightly more smoothed, depending on the regularization scheme used 
(see also Figure [Tot . 



t — 0.5. Figure [O] shows the final result for the magnetic pres- 
sure for the ATHENA run and the three different SPH-MHD imple- 
mentations in GADGET. Visually the results are quite comparable, 
however the GADGET results look slightly more smeared, which 
is the imprint of the underlying SPH and regularization methods. 
This impression is confirmed in figure [141 which shows two cuts 
(y = 0.3125 and y = 0.4277) through the two simulations. Again, 
the black line shows the ATHENA result, the colored lines with the 
error bars are showing the GADGET results. In general there is a 
reasonable agreement, however the SPH-MHD results clearly show 
smoothing of features. The adaptive nature of the SPH-MHD im- 
plementation should allow the central density peak to be resolved 
whereas in ATHENA is can only be resolved by increasing the 



number of grid cells. Never-the-less the SPH-MHD implementa- 
tion seems to converge slower when increasing the resolution (see 
Appendix). 

3.6 General performance of SPH-MHD 

In summary, as shown in the previous sections, an MHD implemen- 
tation in SPH is able to reliable reproduce the results of standard, 
one and two dimensional MHD test problems. We want to stress the 
point that all tests for the SPH-MHD implementation where per- 
formed in a fully three dimensional setup to test the code under re- 
alistic circumstances. Regularization schemes in general are able to 
further suppress the numerical driven growth of div(i3). Although 
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Figure 11. Shown is the resuhing density distribution for the Strong Blast test alt = 0.02. The upper left panel shows the result obtained with ATHENA, 
whereas the upper right panel shows the that from the basic SPH-MHD implementation. The lower left and right panels shows the results obtained with 
the Bsmooth SPH-MHD and the dissipation SPH-MHD implementation respectively. There are no visible difference between the results (see also figure [TH . 
indicating that all SPH-MHD implementations are capable to handle such test situation. 



some optimal numerical values for the regularization schemes can 
be inferred when comparing a suit of different shock tube tests, 
such regularization schemes always introduce small dissipative ef- 
fects, which lead to a slight smearing of sharp features. This has to 
be kept in mind when applying the different SPH-MHD implemen- 
tations to cosmological applications. 



4 COSMOLOGICAL APPLICATION 

The cluster used i n this work is part of a galaxy cluster sample 
jPolag et alJl2008h extracted from a re-simulation of a Lagrangian 
region sele cted from a cosmol ogical, lower resolution DM-only 
simulation jYoshida et al.l200lh . This parent simulation has a box- 



size of 684 Mpc, and assumed a flat ACDM cosmology with 
Q.m = 0.3 for the matter density parameter, Hq = 70 for the Hub- 
ble constant, ftar ~ 0.13 for the baryon fraction and as = 0.9 
for the normalization of the power spectrum. The cluster has a final 
mass of 1.5 X lO"A/0 and was re-simulated at 3 different particle 
masses for the high resolution region. Using the "Zoomed Initial 
Cond itions" (ZIC) technique jXatz & Whitel Il993l ; iTbrmen et al.l 
these regions were re-simulated with higher mass and force 
resolution by populating their Lagrangian volumes with a larger 
number of particles, while appropriately adding additional high- 
frequency modes drawn from the same power spectrum. To opti- 
mize the setup of the initial conditions, the high resolution region 
was sampled with a 16^ grid, where only sub-cells are re-sampled 
at high resolution to allow for quasi abritary shapes of the high 
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Figure 13. The magnetic pressure distribution in the Orszang-Tang Vortex alt = 0.5. The upper left panel shows the result obtained with ATHENA, 

while the upper right panel shows those from the basic SPH-MHD implementation. The lower left and right panels shows the results obtained with the Bsmooth 
SPH-MHD and the dissipation SPH-MHD implementations respectively. Some of the sharp features are more smoothed in the GADGET runs, depending on 
the choice of regularization scheme, but overall the results compare very well (see also figure [T4l 



resolution region. The exact shape of each high-resolution region 
was iterated by repeatedly running dark-matter only simulations, 
until the targeted objects are free of any lower-resolution bound- 
ary particle out to 3-5 virial radii. The initial particle distributions, 
before adding any Zeldov ich displac ement, were taken from a re- 
laxed glass configuration l lWhitelll996 ). The three resolutions used 
correspond to a mass of the dark matter particles of 1.6 x IO^A/q, 
2.5 X 10**Mq and 1.6 x 10** M© for the Ix, 6x and lOx simula- 
tion. The gravitational softening corresponds to 7, 3.9 and 3.2 kpc 
respectively. For simplicity we assumed an initially homogeneous 
magnetic field of 10~^^ G co-moving as also used in previous work 
jPolag et alj[l999ll2002l) . Furthermore we applied the regulariza- 
tion by smoothing the magnetic field in the same way as we did in 



previous work ('Dolagetal."2004','2005') but also tested the effects 
of regularization by artificial dissipation for varying values of q_b. 

Figure [T5] shows a zoom-in from the full cosmological box 
down to the cluster. The structures in the outer parts get less pro- 
nounced due to the decrease in resolution, which is designed to 
capture only the very largest scales of the simulation volume. Each 
panel shows (in clockwise order) a zoom-in by a factor of ten. Fi- 
nally the elongated box in the lower left panel marks the size of 
the observational frame shown on the left. For comparison we pro- 
duced a synthetic Faraday Rotation map from the simulation and 
clipped it to the shape of the actual observations to give an indi- 
cation of the structures resolved by such simulations. We used the 
same linear colorscale for both the observed and the simulated RM 
map, using the highest resoluton simulation {lOx). Note that we 
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Figure 14. Two t = 0.5 cuts throught the pressure in the Orszang-Tang Vortex aly = 0.3125 and at y = 0.4277 (left and right panels, respectively). As 
before, the black line reflect the results obtained with ATHENA and the pink line with the red error bars is obtained with the basic SPH-MHD implementation 
in GADGET. The cuts are choosen for comparison with results from the literature, e.g. lB0rve et alj j20O6l) . 
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Figure 12. Horizontal cut through the Sfrong B/a.vr test (a; = [0.0; 1.0], y = 
0.5) showing the density. The black line is from the ATHENA simulation, 
the pink line with the red error bai's (see figure [Tol reflects the GADGET 
result. The overall behavior is excellent, with only very small differences 
between the two solutions. 

added a constant, galactic foreground signal to the simulated RM 
map to account for the non zero mean in the observed RM map. 
The dynamical range of the simulation spans more than five orders 
of magnitudes in spatial dimension, and the size of the underlying 
box is 6 and 5 time s large r tha n the AMR simulation s presented in 
[Dubois & TevssieJ i tiooj) and lBruggen et al.l l l2005l) , respectively. 
Still the resolution of the underlying dark matter distribution is, re- 
spectively, 2 and 5 times better than these AMR simulations and the 
cluster is resolved with more than one million dark matter particles 
within the virial radius at the lOx resolution. To perform the sim- 
ulation, the lOx resolution run needed ~ 730 CPUh on an AMD 
Opteron cluster. This again is demonstrating the advantages of the 
underlying SPH scheme in making large, cosmological zoomed 
simulations possible. 

To obtain a more quantitative comparison we calculated the 
projected structure function S^^\d) from both the observed and 
synthetic Faraday Rotation maps. 

^'''(Ax, Ay) = {\RM{x,y) - RM{x + Ax,y + Ay)\) ,(37) 

with Ax and Ay being the offsets from a pixel at position 
{x,y). The resulting matrix is then averaged in radial bins d — 



Ax^ + Aj/2 to obtain the structure function. Figure [16] shows 
a comparison of the obtained structure function from the observa- 
tions (black line) and the simulations. For each simulated cluster 
we calculated the synthetic rotation measure maps, clipped accord- 
ingly to the shape of the observed map. To obtain different real- 
izations of the same simulation we produced nine different maps 
where we shifted the clipped region by ±20 kpc in both spacial di- 
rections within the original, cluster centered maps. The thick lines 
mark the mean structure function over these maps, whereas the thin 
lines show the RMS scatter between the different maps. It is clear 
that, due to the additional resolved turbulent velocity field, when 
increasing the resolution the same initial seed fields are amplified 
more. For our lOx simulation the chosen magnetic seed field gets 
amplified to the observed level in our example cluster. Although in- 
creasing the resolution resolve smaller scales in the RM maps, the 
slope of the structure function at small separations, even the lOx 
resolution simulation, gets not as steep as the observed one. This 
confirms the visual impression from the lower left panel in Figure 
[15] which also indicates more structure at small scales in the ob- 
served than the simulated RM map. Pushing the mass resolution by 
one more order of magnitude would probably result in a spacial res- 
olution which should be sufficient to reproduce the observed small 
scale structure in the RM maps, however in such a case we would 
expect to have to start from even smaller seed fields to avoid overes- 
timating the amplitude of the Rotation Measure in the simulations. 
Such a study lies outside of the scope of this paper, as it ultimately 
would lead to questions of the role of magnetic dissipation and vis- 
cosity within the real intra cluster medium. 

In figure [TT] the radial magnetic field profiles are shown for 
the three resolutions comparing the results obtained with the nor- 
mal configuration for cosmological simulations with results where 
we just used the Euler Potential to follow the evolution of the 
magnetic field ign^ring_back reactions. As already noted in ear- 
lier work (Dola g et al.ll2002il . the left panel shows the dependence 
of the amplification of magnetic fields with resolution. In addition, 
the solution obtained with the Euler Potential agrees nicely with 
the B smooth SPH-MHD runs in the outer part of the profiles. It 
is important to note that the increase in amplification of the mag- 
netic field with increasing resolution when using the Euler Poten- 
tial clearly demonstrates that this effect originates from resolving 
more velocity structure, especially driven by the increased amount 
of substructure in the underlying dark matter representation. Thus 
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Figure IS. Zoom into the cluster simulated within the cosmological box. Clockwise, each panel displays a factor 10 increase in imaging magnification, starting 
from the full box (684 Mpc) down to the cluster center (680 kpc). On the very large scale, the density of the dark matter particles are shown, whereas in the 
high resolution region the temperature of the gas is rendered to emphasize the pr esence and dynam ics of the substructure. The last zoom extracts a region of 
the same size of an observed radio jet (3C449) with measured rotation measure jFerettietal.lll999l) . Both, the simulated and the observed map are displayed 
using a linear color-scale based on the minimum and maximum values in the maps. The synthetic RM map is clipped to the shape of the observations. Clearly, 
the simulations are still lacking in resolution, however they do come quite close. 



the result reflects the increased complexity of the structures in the 
density and velocity fields, once the resolution of a cosmological 
simulation is increased. The Euler Potential implementation pro- 
vides an unique possibility to study these effects as they reflect 
the result of integrating the wind-up of the complex flows within 
galaxy clusters, revealing information which can not easily be ob- 
tained in Eulerian schemes. In the central parts, the Bsmooth SPH- 
MHD simulation falls below the solution obtained using the Euler 
Potential. This is easy to understand, because the Euler Potential 
are free from any numerical magnetic dissipation. Additionally the 
magnetic field is strongest in the cluster core and therefore, includ- 
ing the magnetic force in the normal runs will lead to a suppression 
of the amplification. In general the comparison of the two meth- 
ods demonstrate that the amplification of the magnetic field in the 
bsmooth SPH-MHD implementation is not significantly influenced 
by the non-zero div(i3). In addition, although the absolute value 
of the amplification is not converged with resolution, the shape of 
the predicted magnetic field profile appears to be converged. This is 



shown in the right panel of figure[T7] where the profiles are normal- 
ized artificially at large radii to demonstrate the self similar shapes. 
Note that this convergence, as usual for all hydro-dynamical quan- 
tities, is only reached at radii significant larger than the size of the 
gravitational softening, indicated as dashed lines for the lowest (e.g. 
Ix), medium (e.g. 6x) and highest (e.g. lOx) resolution runs. In 
both panels of Figure [T7] we also sho w the results from a cluste r 
simulation using RAMSES, tak en fromlDubois & Tevssieil ( |2008|) . 
and FLASH, taken from Brtig gen et al.l ( 120051) . This comparison 
is over-simplistic, as results are based on simulations of different 
objects and can only be compared with some care. Nevertheless, 
the shape of the radial profiles obtained with RAMSES indicate 
slightly more dissipative effects compared to our Bsmooth SPH- 
MHD implementation, whereas the steeper profile obtained with 
FLASH resembles our results using the Elder Potential implemen- 
tation. 

The situation changes when using artificial magnetic dissipa- 
tion, as shown in figure[T8] The left panel shows the magnetic field 
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Figure 17. The magnetic field profiles for a resolution study of our cluster. Solid lines are obtained with the Bsmooth SPH-MHD implementation, dashed lines 
are for using Euler Potentials. The di fferent line colors indicatin g different res olution. Right panel sh ows the same but normalizing all the profiles in the outer 
part of the cluster. Results taken from lPubois & Tevssiej Ew^) (red line) and lBriiggen et alj i2005l) (diamonds) are also shown. 




Figure 18. The magnetic field profiles obtained for the galaxy cluster using different values of the artificial dissipation. In addition the results using the Euler 
Potentials and the Basic MHD implementation are shown. Right panel shows the same with profiles normalized in the outer part of the cluster. 



profiles for several values of as compared with the profiles for 
the basic SPH-MHD run and that using Euler Potentials. Clearly 
a normal value for artificial magnetic dissipation leads to a large 
dissipation of magnetic field over the simulation time (e.g. close 
to the Hubble time). The right panels show the profiles artificially 
normalized at large radius. Clearly the self similarity of the profiles 
is lost. Therefore it appears that the use of artificially dissipation 
as a regularization scheme is not a good choice for cosmological 
simulations. Additionally it points out that true physical dissipa- 
tion might play an important role in determining the shape of the 
magnetic field profile in galaxy clusters. Here, transport processes, 
cosmic rays, turbulence (especially at unresolved scales) and recon- 
nection of magnetic field lines are not well understood, especially 
within the ICM. As the micro-physical origin of most of them are 
far outside the scales which can be ever reached by cosmological 
simulations, future work will have to include them as approxima- 



tive, sub-grid models, possibly motivated by small scale numerical 
experiments. 



5 CONCLUSIONS 

We presented the implementation of MHD in the cosmological, 
SPH code GADGET. We performed various test problems and dis- 
cussed several instability correction and regularization schemes. 
We also demonstrated the application to cosmological simulations, 
the role of resolution and the role the regularization schemes play 
in cosmological simulations. 
Our main findings are: 

• The combination of many improvements in the SPH im- 
plementatio n, like the correction terms for the variable smooth- 
ing length jSpringel & Hemguisll l2002h as well as the usage of 
the signal velocity in the artificial viscosity ( lMonaghai] [ 19971) to- 
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Figure 19. The magnetic field profiles obtained for the galaxy cluster using different time intervals for the smoothing procedure. In Addition the results for 
using the Euler Potentials and the Basic MHD implementation are shown. Right panel shows the same with profiles normalized in the outer part of the cluster. 
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Figure 16. Shown is the structure function averaged in radial bins (for de- 
tails see text) calculated from the observed (black line) and from the syn- 
thetic Faraday Rotation maps (colored lines). The different colors coiTe- 
spond to different resolutions (Ix, 6x and lOx). The thick lines con'espond 
to the mean calculated over 9 realizations of the maps, whereas the thin lines 
mark the RMS scatter between the different maps (see text for details). 

gether with its generalization to the MHD case jPrice & Monaghai] 
l2004ah improve the handling of magnetic fields in SPH signifi- 
cantly. 

• Correcting the instability by explicitly subtracting the contri- 
bution of a numerical non-zero divergence of the magnetic field 
to the Lorenz for ce from the Maxwell tensor as suggested by 
lB0rve etld] l l200lh seems to perform well. Specifically in three di- 
mensional setups where it seems to work much better than other 
suggestions in the literature. 

• The SPH-MHD implementation performs very well on simple 
shock tube tests as well as on planar test problems. We performed 
all tests in a fully three-dimensional setup and find excellent agree- 
ment of the results obtained with the SPH-MHD implementation 
compared to the results obtained with ATHENA in one or two di- 
mensions. 



• With a convergence study we demonstrate that the SPH-MHD 
results when increasing the resolution are converging to the true 
solution, especially in the sharp features. However, in some regions 
it seems that small but systematic differences converge only very 
slowly to the correct solution. 

• Regularization schemes help to further suppresses noise and 
div(_B) errors in the test simulations, however one has to carefully 
select the numerical parameters to avoid too strong smoothing of 
sharp features. Performing a full set of individual shock tube tests 
allows one to tune the numerical schemes and to determine optimal 
values. However they reflect an optimal choice for problems where 
the local timescales are mostly similar to the global timescale of the 
problem. For cosmological simulations it turns out that regulariza- 
tion by artificial dissipation leads to questionable results, whereas 
the regularization by smoothing the magnetic field (which is ap- 
plied on global timescales) produces reasonable results. 

• The SPH-MHD implementation allows us to perform chal- 
lenging cosmological simulations, covering a large dynamical 
range in length-scales. For galaxy clusters, only the shape of the 
predicted magnetic profiles is, (with the exception of the central 
part of clusters) converged in resolution and in good agreement 
with previous studies. Also the structures obtained in synthetic 
Faraday Rotation maps are in good agreement with previous find- 
ings and compare well with observations. 

The results obtained with artificial dissipation in cosmological 
simulations indicate that physical dissipation could play a crucial 
role in determining the exact shape of the predicted, magnetic field 
profiles in galaxy clusters. Future work, especially when includ- 
ing more physical processes at work in galaxy cluster - as can be 
done easily with our SPH-MHD implementation - will reveal an 
interesting interplay between dynamics of the cluster atmosphere 
and amplification of magnetic fields. Thus having the potential to 
shed light on many, currently unknown aspects of cluster magnetic 
fields, their structure and their evolution. 
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APPENDIX A: CONVERGENCE 

Numerical experiments are normally restricted by the resolution 
one can technically (in terms of computing/memory requirements) 
achieve. Therefore tests as presented in section 3 are usually at 
nominally better resolution than can be obtained in relevant (in 
this case cosmological) simulations. Never-the-less an interesting 
question is, how good do the numerical methods used converge if 
one further increase the resolution? Figure IaTI and lA2l show this 
for Athena and the basic SPH-MHD implementation respectively. 
We repeated the Orszang-Tang Vortex test problem with Athena on 
a 192^, 400^ and 800^ grid. Figure IaT] shows a cut thi'ough the 
density of the Orszang-Tang Vortex, comparing with the result ob- 
tained with the AMR code Ramses ( ITevssiei l2002h . Clearly, the 
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Figure Al. A cut throgh the density for the Orszang-Tang Vortex test (see 
Figure 13/14). Shown in black is the result obtained with Ramses, compared 
to the results obtained with Athena using 3 different resolutions. 
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Figure A2. Same than figure lATI but showing the results obtained with the 
basic SPH-MHD implementation at two different resolutions compared to 
the results obtained with Ramses. 

results obtained with Athena when increasing the resolution ap- 
proaches the results obtained with Ramses. Figure lA2l shows the 
same for setups with 350^ x 5, 700^ x 5 and 1400^ x 5 parti- 
cles. The SPH-MHD implementation also converges towards the 
Ramses results with increasing resolution. However, although the 
central feature is better resolved in the SPH-MHD implementation 
than in the Athena run with comparable resolution, some other fea- 
tures can be seen to converge slower in the SPH-MHD implemen- 
tation when increasing the resolution. Specifically, in some very 
smoothed features there are small but systematic differences be- 
tween the SPH and the true solution. Here the SPH results seem to 
converge only extremely slowly (if at all). 



